# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix Section B
# Appendix Figures 2 & 3 - The relationship between boundaries and crime.

suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# first prepare variables by city
# read in data by city

# Load Atlanta final dataset
load("atl_final.RData")

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv

# Load Austin final dataset
load("aus_final.RData")



# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv

## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv

#### robustness: assessing the relationship between boundaries and crime #### 

# first, evaluating if adjustment attenuates, as opposed to magnifies, the 
# association between the boundary measure and arrests

# atlanta

attm1_atl = lm_robust(lmisdemeanors_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
attm2_atl = lm_robust(lmisdemeanors_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                      data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# austin

attm1_aus = lm_robust(lmisdemeanors_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
attm2_aus = lm_robust(lmisdemeanors_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                      data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# boston

attm1_bos = lm_robust(lfelonies_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
attm2_bos = lm_robust(lfelonies_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
                      data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# chicago

attm1_chi = lm_robust(lmisdemeanors_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

attm2_chi = lm_robust(lmisdemeanors_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                      data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# milwaukee

attm1_mil = lm_robust(lmisdemeanors_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

attm2_mil = lm_robust(lmisdemeanors_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                      data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# seattle

attm1_sea = lm_robust(lmisdemeanors_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                      data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

attm2_sea = lm_robust(lmisdemeanors_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                        hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                      data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

data.frame(
  
  est = c(attm1_atl$coefficients[2],
          attm1_aus$coefficients[2],
          attm1_bos$coefficients[2],
          attm1_chi$coefficients[2],
          attm1_mil$coefficients[2],
          attm1_sea$coefficients[2],
          attm2_atl$coefficients[2],
          attm2_aus$coefficients[2],
          attm2_bos$coefficients[2],
          attm2_chi$coefficients[2],
          attm2_mil$coefficients[2],
          attm2_sea$coefficients[2]),
  
  se = c(attm1_atl$std.error[2],
         attm1_aus$std.error[2],
         attm1_bos$std.error[2],
         attm1_chi$std.error[2],
         attm1_mil$std.error[2],
         attm1_sea$std.error[2],
         attm2_atl$std.error[2],
         attm2_aus$std.error[2],
         attm2_bos$std.error[2],
         attm2_chi$std.error[2],
         attm2_mil$std.error[2],
         attm2_sea$std.error[2]),
  
  pv = c(attm1_atl$p.value[2],
         attm1_aus$p.value[2],
         attm1_bos$p.value[2],
         attm1_chi$p.value[2],
         attm1_mil$p.value[2],
         attm1_sea$p.value[2],
         attm2_atl$p.value[2],
         attm2_aus$p.value[2],
         attm2_bos$p.value[2],
         attm2_chi$p.value[2],
         attm2_mil$p.value[2],
         attm2_sea$p.value[2]),
  
  type = factor(c(rep("No Crime Adjustment", 6),
                  rep("Yes Crime Adjustment", 6)),
                levels = c("No Crime Adjustment",
                           "Yes Crime Adjustment")),
  
  city = factor(rep(c("ATL", "AUS", 'BOS', "CHI", "MIL", "SEA"), 2),
                levels = c("ATL", "AUS", 'BOS', "CHI", "MIL", "SEA"))
  
) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est,
                 col = type),
             position = position_dodge(.4),
             size = 2) + 
  geom_errorbar(aes(x = city,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = type),
                position = position_dodge(.4),
                width = 0,
                size = .4) + 
  geom_errorbar(aes(x = city,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = type),
                position = position_dodge(.4),
                width = 0,
                size = .8) + 
  scale_color_grey(start = .6, end = 0) + 
  labs(x = "City", y = "Boundary Coefficient\n(Standardized)",
       col = "Crime Adjustment?",
       title = 'Adjusting for Crime Always Attenuates the Association Between Racial Boundaries and Arrests') + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 2.5, 
       filename = 'crimeadjust.jpeg')

# ok, now, is crime downstream of boundaries? 

# atlanta

ds1 = lm_robust(scale(lpropertycrime) ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
ds2 = lm_robust(scale(lviolentcrime) ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# austin

ds3 = lm_robust(scale(lpropertycrime) ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
ds4 = lm_robust(scale(lviolentcrime) ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# boston

ds5 = lm_robust(scale(lcrime) ~ bos_white_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# chicago

ds6 = lm_robust(scale(lpropertycrime) ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

ds7 = lm_robust(scale(lviolentcrime) ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# milwaukee

ds8 = lm_robust(scale(lpropertycrime) ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

ds9 = lm_robust(scale(lviolentcrime) ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                  hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# seattle

ds10 = lm_robust(scale(lpropertycrime) ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                 data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

ds11 = lm_robust(scale(lviolentcrime) ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol, 
                 data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

data.frame(
  
  est = c(ds1$coefficients[2],
          ds2$coefficients[2],
          ds3$coefficients[2],
          ds4$coefficients[2],
          ds5$coefficients[2],
          ds6$coefficients[2],
          ds7$coefficients[2],
          ds8$coefficients[2],
          ds9$coefficients[2],
          ds10$coefficients[2],
          ds11$coefficients[2]),
  
  se = c(ds1$std.error[2],
         ds2$std.error[2],
         ds3$std.error[2],
         ds4$std.error[2],
         ds5$std.error[2],
         ds6$std.error[2],
         ds7$std.error[2],
         ds8$std.error[2],
         ds9$std.error[2],
         ds10$std.error[2],
         ds11$std.error[2]),
  
  pv = c(ds1$p.value[2],
         ds2$p.value[2],
         ds3$p.value[2],
         ds4$p.value[2],
         ds5$p.value[2],
         ds6$p.value[2],
         ds7$p.value[2],
         ds8$p.value[2],
         ds9$p.value[2],
         ds10$p.value[2],
         ds11$p.value[2]),
  
  outcome = c("Property", "Violent", "Property", "Violent", "All", "Property", "Violent",
              "Property", "Violent", "Property", "Violent"),
  
  city = c("ATL", "ATL", "AUS", "AUS", 'BOS', "CHI", "CHI", "MIL", "MIL", "SEA", "SEA")
  
) %>% 
  ggplot() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  geom_point(aes(x = city, y = est,
                 col = outcome),
             position = position_dodge(.4),
             size = 2) + 
  geom_errorbar(aes(x = city,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = outcome),
                position = position_dodge(.4),
                width = 0,
                size = .4) + 
  geom_errorbar(aes(x = city,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = outcome),
                position = position_dodge(.4),
                width = 0,
                size = .8) + 
  scale_color_grey(start = .6, end = 0) + 
  labs(x = "City", y = "Boundary Coefficient\n(Standardized)",
       col = "Outcome",
       title = 'Racial Boundaries are Positively Associated With Crime') + 
  theme_tufte()

ggsave(plot = last_plot(), width = 8, height = 2.5, 
       filename = 'rbcrimerel.jpeg')
